Use R package UScensus2010county to complete the following tasks: (20 pt.)
Plot a map of New York counties using the plot function.
if(!require('UScensus2010')) install.packages('UScensus2010')
library(UScensus2010)
if(!require('UScensus2010county')) install.county('osx')
library(UScensus2010county)
data('new_york.county10')
shp <- new_york.county10
plot(shp)
Plot a map of New York counties using the qtm function.
if(!require('tmap')) install.packages('tmap')
library(tmap)
qtm(shp,fill = 'P0010001')
How many counties in New York State?
df <- shp@data
df
nrow(df)
[1] 62
What’s the 3-digit fips code of Broome County?
str(df)
'data.frame': 62 obs. of 462 variables:
$ NAME10 : chr "Broome" "Seneca" "Schenectady" "Niagara" ...
$ NAMELSAD10: chr "Broome County" "Seneca County" "Schenectady County" "Niagara County" ...
$ state : chr "36" "36" "36" "36" ...
$ county : chr "007" "099" "093" "063" ...
$ fips : chr "36007" "36099" "36093" "36063" ...
$ P0010001 : num 200600 35251 154727 216469 159429 ...
$ P0020001 : num 200600 35251 154727 216469 159429 ...
$ P0020002 : num 0 0 0 0 0 0 0 0 0 0 ...
$ P0020003 : num 0 0 0 0 0 0 0 0 0 0 ...
$ P0020004 : num 0 0 0 0 0 0 0 0 0 0 ...
$ P0020005 : num 0 0 0 0 0 0 0 0 0 0 ...
$ P0020006 : num 200600 35251 154727 216469 159429 ...
$ P0030001 : num 200600 35251 154727 216469 159429 ...
$ P0030002 : num 176444 32591 123211 191673 139529 ...
$ P0030003 : num 9614 1607 14710 14851 10338 ...
$ P0030004 : num 396 104 575 2285 385 ...
$ P0030005 : num 7065 244 4960 1823 3517 ...
$ P0030006 : num 82 2 105 62 34 34 213 16 8 19 ...
$ P0030007 : num 1912 240 5364 1108 1783 ...
$ P0030008 : num 5087 463 5802 4667 3843 ...
$ P0040001 : num 200600 35251 154727 216469 159429 ...
$ P0040002 : num 193822 34299 145900 211775 153349 ...
$ P0040003 : num 6778 952 8827 4694 6080 ...
$ P0050001 : num 200600 35251 154727 216469 159429 ...
$ P0050002 : num 193822 34299 145900 211775 153349 ...
$ P0050003 : num 173074 31999 119409 188907 136555 ...
$ P0050004 : num 8850 1513 13528 14511 9592 ...
$ P0050005 : num 328 96 445 2135 283 ...
$ P0050006 : num 7019 238 4917 1807 3469 ...
$ P0050007 : num 60 2 81 55 25 29 137 14 5 15 ...
$ P0050008 : num 242 47 2708 185 241 ...
$ P0050009 : num 4249 404 4812 4175 3184 ...
$ P0050010 : num 6778 952 8827 4694 6080 ...
$ P0050011 : num 3370 592 3802 2766 2974 ...
$ P0050012 : num 764 94 1182 340 746 ...
$ P0050013 : num 68 8 130 150 102 183 1000 55 27 20 ...
$ P0050014 : num 46 6 43 16 48 46 467 14 7 8 ...
$ P0050015 : num 22 0 24 7 9 5 76 2 3 4 ...
$ P0050016 : num 1670 193 2656 923 1542 ...
$ P0050017 : num 838 59 990 492 659 ...
$ P0060001 : num 206107 35732 161052 221479 163567 ...
$ P0060002 : num 181009 33034 127516 195911 143012 ...
$ P0060003 : num 12499 1812 17992 17632 12620 ...
$ P0060004 : num 1664 281 1769 3859 1273 ...
$ P0060005 : num 8193 291 6587 2407 4293 ...
$ P0060006 : num 243 8 391 158 123 165 683 65 34 58 ...
$ P0060007 : num 2499 306 6797 1512 2246 ...
$ P0070001 : num 206107 35732 161052 221479 163567 ...
$ P0070002 : num 198382 34718 151125 216232 156746 ...
$ P0070003 : num 176969 32386 123004 192737 139499 ...
$ P0070004 : num 11278 1705 16273 17058 11487 ...
$ P0070005 : num 1445 269 1465 3589 1058 ...
$ P0070006 : num 8072 285 6436 2354 4205 ...
$ P0070007 : num 191 8 327 141 100 122 439 51 25 42 ...
$ P0070008 : num 427 65 3620 353 397 ...
$ P0070009 : num 7725 1014 9927 5247 6821 ...
$ P0070010 : num 4040 648 4512 3174 3513 ...
$ P0070011 : num 1221 107 1719 574 1133 ...
$ P0070012 : num 219 12 304 270 215 ...
$ P0070013 : num 121 6 151 53 88 ...
$ P0070014 : num 52 0 64 17 23 43 244 14 9 16 ...
$ P0070015 : num 2072 241 3177 1159 1849 ...
$ P0080001 : num 200600 35251 154727 216469 159429 ...
$ P0080002 : num 195513 34788 148925 211802 155586 ...
$ P0080003 : num 176444 32591 123211 191673 139529 ...
$ P0080004 : num 9614 1607 14710 14851 10338 ...
$ P0080005 : num 396 104 575 2285 385 ...
$ P0080006 : num 7065 244 4960 1823 3517 ...
$ P0080007 : num 82 2 105 62 34 34 213 16 8 19 ...
$ P0080008 : num 1912 240 5364 1108 1783 ...
$ P0080009 : num 5087 463 5802 4667 3843 ...
$ P0080010 : num 4703 445 5296 4347 3577 ...
$ P0080011 : num 2181 181 2242 2165 1817 ...
$ P0080012 : num 843 160 467 1094 588 ...
$ P0080013 : num 797 36 847 397 576 ...
$ P0080014 : num 55 2 16 45 27 36 134 19 8 13 ...
$ P0080015 : num 318 47 348 230 226 ...
$ P0080016 : num 147 5 198 202 91 175 523 31 12 19 ...
$ P0080017 : num 68 0 111 42 44 65 348 10 10 9 ...
$ P0080018 : num 29 0 64 16 8 13 53 2 4 5 ...
$ P0080019 : num 120 3 271 77 93 160 953 21 7 15 ...
$ P0080020 : num 18 1 72 22 9 39 199 15 2 3 ...
$ P0080021 : num 2 0 3 2 3 8 12 4 1 1 ...
$ P0080022 : num 8 2 107 14 13 41 151 3 1 2 ...
$ P0080023 : num 26 1 28 10 11 15 70 8 4 9 ...
$ P0080024 : num 80 7 420 29 62 51 620 9 4 7 ...
$ P0080025 : num 11 0 102 2 9 18 95 5 1 2 ...
$ P0080026 : num 354 18 489 298 238 360 819 54 46 82 ...
$ P0080027 : num 202 7 236 191 135 225 380 33 24 50 ...
$ P0080028 : num 67 1 36 30 17 33 66 1 4 7 ...
$ P0080029 : num 5 2 3 2 7 3 10 2 4 0 ...
$ P0080030 : num 23 5 56 24 29 36 87 3 2 3 ...
$ P0080031 : num 12 0 7 16 18 7 29 2 6 4 ...
$ P0080032 : num 0 0 4 1 2 2 7 0 1 0 ...
$ P0080033 : num 5 1 3 5 4 11 20 4 1 2 ...
$ P0080034 : num 16 0 11 12 3 17 20 3 1 7 ...
$ P0080035 : num 9 1 7 2 4 8 91 1 2 4 ...
$ P0080036 : num 2 0 5 2 2 2 21 1 1 0 ...
$ P0080037 : num 2 0 8 5 1 7 20 0 0 3 ...
[list output truncated]
df[df$NAME10 == 'Broome',]['fips']
Compute descriptive statistics of the population column (P0010001), including total, minimum, maximum, mean, median, and skewness.
if(!require(moments)) install.packages('moments')
library(moments)
population <- df$P0010001
total <- sum(population)
minimum <- min(population)
maximun <- max(population)
mean <- mean(population)
median <- median(population)
skewness <- skewness(population)
total
[1] 19378102
minimum
[1] 4836
maximun
[1] 2504700
mean
[1] 312550
median
[1] 91301
skewness
[1] 2.579801
Create a histogram and a boxplot of the population.
hist(population)
boxplot(population)
Use R package UScensus2010tract to complete the following tasks: (20 pt.)
if(!require('UScensus2010tract')) install.tract('osx')
library(UScensus2010tract)
Plot a map of New York census tracts using the plot function.
data(new_york.tract10)
shp <- new_york.tract10
plot(shp)
Compute the total population based on census tracts.
df <- shp@data
sum(df$P0010001)
[1] 19378102
df
Select all census tracts in Broome County and plot the map.
plot(shp[df$county == '007',])
What’s the total population of Broome County?
broome <- shp[df$county == '007',]
broomedf <- broome@data
sum(broomedf$P0010001)
[1] 200600
Create a histogram and a boxplot of population based on census tracts of Broome County.
hist(broomedf$P0010001)
boxplot(broomedf$P0010001)
Select the first five columns of the shapefile of Broome County census tract; add a new population ratio column (= census tract population / county population); save the new shapefile to the hard drive.
broome_5 <- broome[,1:5]
df <- broome_5@data
df$popratio <- df$P0010001 / sum(broomedf$P0010001)
broome_5@data <- df
if(!require(rgdal)) install.packages('rgdal')
library(rgdal)
getwd()
[1] "C:/Users/hwang164/Documents/GEOG-533/Labs/Lab_09"
writeOGR(broome_5,dsn = '.', layer = "broome_5", driver = 'ESRI Shapefile')
Error in writeOGR(broome_5, dsn = ".", layer = "broome_5", driver = "ESRI Shapefile") :
layer exists, use a new layer name
Use R packages raster and ncdf4 to complete the following tasks: (20 pt.)
if(!require('raster')) install.packages('raster')
library(raster)
if(!require('ncdf4')) install.packages('ncdf4')
Loading required package: ncdf4
there is no package called 㤼㸱ncdf4㤼㸲Installing package into 㤼㸱C:/Users/hwang164/Documents/R/win-library/3.5㤼㸲
(as 㤼㸱lib㤼㸲 is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/3.5/ncdf4_1.16.1.zip'
Content type 'application/zip' length 2598248 bytes (2.5 MB)
downloaded 2.5 MB
package ‘ncdf4’ successfully unpacked and MD5 sums checked
The downloaded binary packages are in
C:\Users\hwang164\AppData\Local\Temp\RtmpOegO0P\downloaded_packages
library(ncdf4)
package 㤼㸱ncdf4㤼㸲 was built under R version 3.5.3
Load the multi-band raster dataset “NDVI.nc” into R.
ndvi.rb <- brick('NDVI.nc')
Get the basic information about the dataset, including the number of rows, columns, cells, and bands; spatial resolution, extent, bounding box, and projection.
nrow(ndvi.rb)
[1] 360
ncol(ndvi.rb)
[1] 720
ncell(ndvi.rb)
[1] 259200
nlayers(ndvi.rb)
[1] 36
res(ndvi.rb)
[1] 0.5 0.5
extent(ndvi.rb)
class : Extent
xmin : -180
xmax : 180
ymin : -90
ymax : 90
bbox(ndvi.rb)
min max
s1 -180 180
s2 -90 90
projection(ndvi.rb)
[1] "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
Aggregate all bands to generate a mean NDVI raster and a maximum NDVI raster; save the two new raster datasets to the hard drive.
mean <- mean(ndvi.rb)
max <- max(ndvi.rb)
plot(mean)
plot(max)
writeRaster(mean, filename = 'mean_ndvi')
class : RasterLayer
dimensions : 360, 720, 259200 (nrow, ncol, ncell)
resolution : 0.5, 0.5 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : C:/Users/hwang164/Documents/GEOG-533/Labs/Lab_09/mean_ndvi.grd
names : layer
values : 0, 0.8537528 (min, max)
writeRaster(max, filename = 'max_ndvi')
class : RasterLayer
dimensions : 360, 720, 259200 (nrow, ncol, ncell)
resolution : 0.5, 0.5 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : C:/Users/hwang164/Documents/GEOG-533/Labs/Lab_09/max_ndvi.grd
names : layer
values : 0, 0.9681 (min, max)
Plot the maps of monthly NDVI of the year 2001
Create histograms of monthly NDVI of the year 2001.
hist(ndvi2001)
size changed to n because it cannot be larger than n when replace is FALSE39% of the raster cells were used (of which 42% were NA). 58093 values used.
Plot the NDVI map of July 2000; click any location with data on the map and retrieve the NDVI time series for all years; plot the NDVI time series of the selected location.
values <- click(ndvi.rb, n = 1, xy = FALSE)
Error in graphics::locator(1, type, ...) :
plot.new has not been called yet